arXiv: 1509.0307 lv2 [physics.chem-ph] 5 May 2016 


Thermal density functional theory: Time-dependent linear response and approximate 
functionals from the fluctuation-dissipation theorem 

Aurora Pribram-Jones,^A Paul E. Grabowski,^ and Kieron Burke^’^ 

^Lawrence Livermore National Laboratory, Livermore, CA 94550 
^Department of Chemistry, University of California, Berkeley, CA 94720 
^Department of Physics and Astronomy, University of California, Irvine, CA 92697 
Department of Chemistry, University of California, Irvine, CA 92697 
(Dated: May 9, 2016) 

The van Leeuwen proof of linear-response time-dependent density functional theory (TDDFT) is 
generalized to thermal ensembles. This allows generalization to finite temperatures of the Gross- 
Kohn relation, the exchange-correlation kernel of TDDFT, and fluctuation dissipation theorem 
for DFT. This produces a natural method for generating new thermal exchange-correlation (XC) 
approximations. 
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Kohn-Sham density functional theory (KS-DFT) is 
a popular and well-established approach to electronic 
structure problems in many areas, especially materi¬ 
als science and chemistry [1]. The Kohn-Sham method 
imagines a fictitious system of non-interacting fermions 
with the same density as the real system[2] and from 
which the ground-state energy can be extracted. Only 
a small fraction of the total energy, called the exchange- 
correlation (XC) energy, need be approximated to solve 
any ground-state electronic problem[l], and modern ap¬ 
proximations usually produce sufficient accuracy to be 
useful [3]. The advent of TDDFT generalized this method 
to time-dependent problems [4]. Limiting TDDFT to 
linear-response yields a method for extracting electronic 
excitations [5, 6], once another functional, the XC kernel, 
is also approximated. 

But there is growing interest in systems in which the 
electrons are not close to zero temperature. Warm dense 
matter (WDM) is partially ionized, solid-density mat¬ 
ter having a temperature near the Fermi energy. It has 
wide-ranging applications including the astrophysics of 
giant planets and white dwarf atmospheres [7-14], cheap 
and ultra-compact particle accelerators and radiation 
sources[15-17], and the eventual production of clean, 
abundant energy via inertial confinement fusion[18, 19]. 
One of the most successful methods for simulating equi¬ 
librium warm dense matter combines DFT[2, 20] and 
molecular dynamics[21] to capture quantum mechani¬ 
cal effects of WDM electrons and the classical behavior 
of ions [7-14, 22-24]. Such simulations use the Mermin 
theorem[25] to generate a KS scheme at finite tempera¬ 
ture, defined to generate the equilibrium density and free 
energy. In practice, the XC free energy is almost always 
approximated with a ground-state approximation, but 
formulas for thermal corrections are being developed [26- 
30]. 

Many processes of interest involve perturbing an equi¬ 
librium system with some time-dependent (TD) pertur¬ 
bation, such as a laser field[31] or a rapidly moving 


nucleus as in stopping power [32-34]. Of great interest 
within the WDM community are calculations of spectra, 
dynamic structure factors, and the flow of energy be¬ 
tween electrons and ions [35-38]. Spectra expose a mate¬ 
rial’s response to excitation by electromagnetic radiation, 
which would facilitate experimental design and analysis. 
Dynamic structure factors can be related to the x-ray 
scattering response, which is being developed as a tem¬ 
perature and structural diagnostic tool for WDM[39]. 
Thus it would appear that a TD version of the Mer¬ 
min formalism is required. A theorem is proven in Li 
et al.[40, 41], but the formalism assumes the temper¬ 
ature is fixed throughout the process, and so cannot 
describe e.g., equilibration between electrons and ions. 
Moreover, the proof requires the Taylor expansion of the 
perturbing potential as a function of time, just as in the 
Runge-Gross (RG) theorem[4]. This can be problematic 
for initial states with cusps [42], such as at the nuclear 
centers. (Recent efforts[43, 44] have focused on avoiding 
these complications at zero temperature.) Finally, the 
RG proof requires invocation of a boundary condition to 
complete the one-to-one correspondence between density 
and potential[45], which create subtleties when applied 
to extended systems [46]. 

In the present work, we prove the RG theorem at finite 
temperature within linear response by generalizing the el¬ 
egant linear response proof of van Leeuwen [43] to thermal 
ensembles. Our proof avoids several of the drawbacks 
mentioned above, while still providing a solid ground¬ 
ing to much of WDM theoretical work. We then de¬ 
fine the exchange-correlation kernel at finite temperature 
and generalize the Gross-Kohn equation. Finally, we ex¬ 
tend the fluctuation-dissipation theorem of ground-state 
DFT to finite temperatures, and show how this provides 
a route to equilibrium free energy XC approximations. 

Consider a system of electrons in thermal and parti¬ 
cle equilibrium with a bath at some temperature, r, and 
with static equilibrium density '^’’’(r). The system ex¬ 
tends throughout space with a finite average density, i.e.. 
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the thermodynamic limit has been taken. The limit of 
isolated atoms or molecules is achieved by then taking 
the separation between certain nuclei to infinity. In this 
sense, no surface boundary condition need be invoked [45], 
as the density never quite vanishes, while the average 
particle number per atom or molecule molecule is finite. 
These electrons are perturbed at t = 0 by a potential 
Sv{r,t) that is Laplace-transformable. To avoid complex 
questions of equilibration, we consider only the linear re¬ 
sponse of the system, so that the perturbation does not 
affect the temperature of the system as, e.g.. Joule heat¬ 
ing is a higher order effect [47]. The Kubo response for¬ 
mula for the density change in response to Sv is 

(5n^(r,s)= J d^r' x'"{r,r',s)5v{r',s), (1) 


where the Laplace transform 

Sv{r,s)= / Sv{r,t) (2) 

Jo 

is assumed to exist for all s > 0. Within the grand canon¬ 
ical ensemble[48, 49], the equilibrium density-density re¬ 
sponse function is [50]: 


X 


■(r,r',s) = Wi 


An[*(r)An[-(r') 


S — lUJi 


-bc.c., (3) 


where 


Auy W = (*|^(r)|j) - Sij n^(r) (4) 

are matrix elements of the density fluctuation operator. 
The energy-ordered indices i,j run over all many-body 
states (both bound and continuum[51]) with all parti¬ 
cle numbers, but An[j vanishes unless W = Nj. The 
transition frequencies ujji = Ej — Ei, and the statistical 
weights Wi are thermal occupations for the equilibrium 
statistical operator r”" = (^*1 obey < Wj 

if Ei > Ej and W = Nj. This condition is satisfied by 
the grand canonical ensemble of common interest with 
Wi = exp[-(£'i - / t]/ - iJLNi)/T]. 

We also need the (Laplace-transformed) one-body po¬ 
tential operator: 

i5I4(s) = J d^rh{r)Sv{r,s), (5) 

and its matrix elements: 

SV.,{s) = {2\6ns)\j)- ( 6 ) 

Its expectation value is 

6V^{s) ='^WiVu{s) = J dJrTE(r)5v{v,s), (7) 

i 

so that matrix elements of its fluctuations are 


AV^^{s)=6V,,{s)-6,,SV^s). ( 8 ) 


Then consider the expectation value: 

mJ'{s) = J d^rSn'^{r,s)6v{r,s). (9) 

Inserting Eq. (1) and using the definitions, we find 


-p ■ 


This is rearranged as 


^(s) = -2 


oo oo 
i=0 j=i+l ^ ^ ^3^ 


I AW(s) |2 


( 11 ) 


We have ordered all states by energy regardless of parti¬ 
cle number here for simplicity, though this is not strictly 
necessary since different particle number subsystems do 
not interact. For now, we assume no degeneracies. 
Then the above expression, m*'(s), vanishes only if ev¬ 
ery AV^j{s) does for i ^ j because of our assumption 
that (wi — Vj)u!ji > 0 a i ^ j. 

The usual statement of the RG theorem is that no two 
potentials that differ by more than an inconsequential 
function of time alone can give rise to the same density 
(for fixed statistics, interparticle interaction, and initial 
state[4]). Imagine two such perturbations exist, yielding 
the same density response. Since, in linear response, the 
density response is proportional to the perturbation, we 
can subtract one from the other, and the statement to be 
proved is that there is no non-trivial perturbation with 
zero density response. If it did exist, then rrE{s) would 
vanish and our algebra shows that every AVjj{s) with 
i ^ j would also. Finally, 

Ni 

^6v{rk,s)^j{ri ... tatJ = X 

k—1 i 

. 

which can be proven by integrating over all coordinates 
with 4*^. Then, as AV^j{s) = SVij{s) for i ^ j, and must 
vanish if there is no density response, the sum on the 
right of Eq. (12) collapses to just the j-th term, showing 
that 5v{r, s) must be spatially independent. 

We can also include a finite number (M) of degener¬ 
ate excited eigenstates. (For the complications involved 
when the ground-state is degenerate, see Ref. [52]). For 
such states, oJij = 0 and the argument above no longer 
implies SVij(s) vanishes, as the perturbation couples de¬ 
generate states within the same subspace. But simply 
choose at least M points in the 3A-dimensional coordi¬ 
nate space that are not on any nodal hypersurface of the 
degenerate subspace. Then the only solution to Eq. (12) 
is again that Jn(r, s) must be independent of r. 

Thus we have generalized the van Leeuwen proof to 
thermal ensembles, even with finite degeneracies among 
excited states. Our proof applies to any ensemble with 
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weights that monotonically decrease with increasing en¬ 
ergy for each particle number[53, 54]. This avoids compli¬ 
cations caused by cusps in initial wavefunctions[42, 55]. 
Extension to spatially periodic potentials is straightfor¬ 
ward, as no boundary condition [45] was invoked [46]. 

In order for the above result to be of practical use, 
we consider the KS scheme for finite-temperature, time- 
dependent systems and provide a method for generating 
XC approximations. We assume the equilbrium Mermin- 
Kohn-Sham (MKS)[2, 25] potential exists. At this point, 
we switch to using the more familiar Fourier-transform 
notation, but in fact all results and definitions apply 
only to Laplace-transformable perturbations. (In prac¬ 
tice, this distinction rarely matters, but occasional formal 
difficulties arise if this restriction is not made, see Ref. 
[56] and Sec. 3.2 of Ref. [43].) First we generalize the 
Gross-Kohn response formula[57] to thermal ensembles. 
Define 


A simple approximation is then the thermal adiabatic lo¬ 
cal density approximation (thALDA), in which the ther¬ 
mal XC kernel is approximated using the XC free energy 
density per particle for a finite-temperature uniform gas. 


thALDA 
J XC 


[n](r,r',a;) 


72 T.unif / \ 

Q Qxc {nj 
<Pn 



r'), (19) 


which ignores its nonlocality in space and time, and could 
be used to generalize ALDA calculations of excitations in 
metals and their surfaces[61]. 

We next deduce the fluctuation-dissipation theorem for 
MKS thermal DFT calculations. This allows us to con¬ 
nect the response function and the Coulomb interaction 
through the dynamical structure factor [62]. In the MKS 
scheme, the XC contributions to the free energy are de¬ 
fined via 


X^(r,r',a;) = 

y- ^ f (j|n(r)|fc)(fc|n(r')|j) _ (j|n(F)|fc)(fc|h(r)|j) 
I w - UJkj +iT] to + ujkj + ir] 

jk 

(13) 

where 77 —>• 0+[58]. 

Because of our proof of one-to-one correspondence, we 
can invert the response function (excluding a constant), 
and write 


(X")-'( 12 ) 


du(l) 
(5n(2) ’ 


(14) 


where 1 denotes the coordinates r, t, and 2 another 
pair[59]. The standard definition of XC is: 


7;s(l) = u(l)-|-U h( 1)+Uxc(l), (15) 


where v^, is the one-body potential of the non-interacting 
KS system and Wh is the Hartree potential [60]. Differen¬ 
tiating with respect to n( 2 ), this yields 

{xir" (12) = (x")”' (12) + /h(12) + /^c(12), (16) 


which defines the XC kernel at finite temperature, where 
Xl is the KS response function [58] and the traditionally 
defined Hartree contribution is simply 


/h(12) 


5{tl — t2) 

I ri - r 2 I ■ 


(17) 


This follows the definition within the Mermin 
formalism[25] (but see Refs. [48] and [53] for alter¬ 
native choices and their consequences). Inverting yields 
the thermal Gross-Kohn equation[57]: 


X"(12)=x^(12)-b J d3d4x;(13)/^xc(34)x"(42). (18) 


A'^[n] = T[n] + t4e[n] -|- V[n] — T5'[n] (20) 

= Ts W + U[n] + V[n] - r5s[n] -|- (21) 

By subtraction. 

Ale [«] = ^xc [n] + Ule [?^] - tS'xc [n] (22) 


where T denotes kinetic, U denotes potential, and S 
entropic components. Using many-body theory, the 
density-density response function determines the poten¬ 
tial contribution to correlation [63, 64], just as in the 
ground state [65]: 


ui = 


Ke[f"[n]] - UeeKN] 





(23) 

duj . 1 

— coth 

^ w A AAx^(r,r',a;) 

27r ' 

^2tJ |r-r'| 


(24) 


where Ax’" = x’" ~ Xs- introducing a coupling- 

constant A while keeping the density fixed, the thermal 
connection formula[ 66 ] yields 

Alin] = 

SAx'"'K](r,r',a;) 
|r - r'l 
(25) 

where the scaled density is n^{r) = x^n(xv) and 7 = 
yj t' jr. This is exact, but only if the exact thermal 
XC kernel is used, as defined by Eq. (16). If the ker¬ 
nel is omitted, the result is the thermal random-phase 
approximation[67]. 

Next, we discuss the many applications of Eq. (25). 

There has been tremendous progress in implementing 
and testing the random phase approximation for calcu¬ 
lating the XC energy in ground-state calculations and 


lim - 

T"—yOO 2 


^4 fdr 


j f dio I ^ 
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27r 
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such calculations, while more expensive than standard 
DFT, are becoming routine [68-70]. Our results pro¬ 
vide a thermal generalization that could likewise be used 
to generate new thermal XC approximations for equi¬ 
librium WDM calculations. At finite temperature, the 
XC hole fails to satisfy the simple sum rules [71] that 
have proven so powerful in constructing ground-state 
approximations [72]. But our formula uses instead the 
XC kernel. Inserting Eqs (18,19) into Eq. (25) yields 
thALDA-RPA, a new approximation to the equilibrium 
correlation energy, that can be applied to any system. 
Another, simpler approximation is ALDA, in which only 
the zero-temperature XC energy is used in the kernel. 
Both can be relatively easily evaluated for a uniform gas, 
and the resulting a^^{rs.) found from Eq. (25) compared 
with an accurate parametrization[27]. Even in the uni¬ 
form gas, thALDA is an approximation because both the 
q- and w-dependence of the true /Jc are missing; thus 
the efficacy of these approximations can be tested on the 
uniform case. 

Next we discuss which known exact conditions on the 
zero-temperature kernel apply to the thermal kernel, and 
which do not. Because the equilibrium solution is a mini¬ 
mum of the thermal free-energy functional, the zero-force 
theorem [64] 

J Sr J ST'Sir)n^r')f^^{r,r',Lo) = 0 (26) 

should be satisfied and the kernel should be symmetric in 
its spatial arguments. However, any simple formula for a 
one-electron system[73] is not true at finite temperature, 
as the particle number is only an average in the grand 
canonical ensemble[49, 71]. 

A last set of conditions is found by considering the 
coupling-constant dependence in DFT. A parameter A is 
introduced that multiplies the electron-electron interac¬ 
tion, while keeping the density constant. Because of sim¬ 
ple scaling relations, the A-dependence can be shown to 
be determined entirely by coordinate scaling of the den¬ 
sity as in Eq. (25), i.e., determined by the functional it¬ 
self, evaluated at different densities. This is used in both 
ground-state DFT[74] and in time-dependent DFT[75], 
and has been generalized to the thermal case[66, 76]. 
Although the thermal connection formula does not re¬ 
quire this relation for the response function, it is useful 
in many contexts. From the Lehmann representation[50] 
of x’’ [63] , we find the A-dependent response function sat¬ 
isfies: 

X’'’''[n](r,r',w) = A'‘[’^i/A](Ar, Ar', w/A^). (27) 

Insertion into the definition of /xc yields: 

= AV^c^ K/A](Ar, Ar',aj/A^), (28) 

and the potential perturbation scales as: 

[«](!■,w) = A^^uxc^ K/A](Ar,w/A2). (29) 


Insertion of the scaling relation for the kernel into the 
thermal connection formula yields a more familiar analog 
to the ground-state formula. 

The exchange kernel must scale linearly with coupling 
constant, so Eq. (28) produces a rule for scaling of the 
exchange kernel: 

/][;[n^](r,r',w) = 7 /x^^ [n]( 7 r,yr',w/y^). (30) 

Because the poles in /xc are A-dependent, we expect 
pathologies similar to those in zero-temperature TDDFT 
if the exact frequency-dependent // is used in Eq. 
(25) [77]. But adiabatic EXX (AEXX), not including 
frequency-dependence, produces a well-defined approx¬ 
imation to the thermal free energy in which the kernel is 
non-local. This and the other proposed approximations 
above could prove useful in WDM simulations when ther¬ 
mal XC effects are relevant (but see [78] for discussion of 
the subtleties involved in thermal XC approximations). 

In conclusion, we have generalized the proofs and con¬ 
structions of TDDFT within the linear response formal¬ 
ism to thermal ensembles, including those containing a 
finite number of degeneracies. We have avoided ambigu¬ 
ities about the relative perturbative and thermal equili¬ 
bration time scales, allowed for degenerate excited states 
more common in finite-temperature ensembles, avoided 
invoking boundary conditions and the requirement of 
Taylor expandability, and provided firm footing for finite- 
temperature, time-dependent KS-DFT in the linear re¬ 
sponse regime. Definition of relevant KS quantities led 
to description of their properties under scaling. Further, 
we have shown that these quantities, in combination with 
the thermal connection formula, produce new routes to 
thermal DFT approximations for use in equilibrium MKS 
calculations. Implementation and tests of these approxi¬ 
mations is ongoing. 
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